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1. Introduction 



Magnetic fields play an important role in the behavior of plasmas and are thought mediate 
important effects like dynamos in the core of planets and the formation of jets in active galactic 
^ I nuclei and gamma ray bursts; induce a variety of magnetic instabilities; realize solar flares, etc. 

(see e.g. p^ ). Understanding the role of magnetic fields in these and other phenomena have 
spurred through the years many efforts to obtain solutions of the magnetohydrodynamic equa- 
tions. The non-linear nature of these equations limits the understanding that can be gained 
O I in a particular problem via analytical techniques. This implies that solutions for complex sys- 

■ tems must be obtained by numerical means and a suitable numerical implementation must be 
^ . constructed for this purpose. Such implementation must be able to evolve the solution to the 

future of some initial configuration and guarantee its quality. A delicate, subsidiary quantity, 
O ' can be monitored in part to estimate this. This quantity is the "monopole constraint" diB^ 
which must be zero at the analytical level for a consistent solution. This quantity is not a part 
of the main variables, rather it is a derived quantity which should be satisfied by a true physi- 
cal solution. In practice, unless a numerical implementation of the MHD equations is carefully 
^ I designed this constraint can be severely violated. This, in turn, signals (and is sometimes the 

■ cause) of a degrading numerical solution. For these reasons, several approaches have been in- 
vestigated and developed for guaranteeing a controlled behavior of this quantity. One such 
approach is known as the constraint transport technique which adopts a particular algo- 
rithm that staggers the variables appropriately to ensure the satisfaction of the constraint at 
round-off level within Finite Difference and Finite Elements techniques. This approach has 
been quite successful in a number of applications across different disciplines and particularly 
relevant in astrophysics applications [5|l6|l7|l8|9|10|ll|12j . However, by design it imposes limits 
on the algorithmic options available to an implementation. This fact can be at odds, or intro- 
duce complications, with applications where adaptive mesh refinement is required and/or ad- 
vanced numerical techniques (that exploit useful properties of the equations) are adopted. An 
alternative approach, which controls the constraint at truncation-error level maintains com- 
plete freedom in the numerical techniques to be adopted. This approach, referred to as diver- 
gence cleaning puts the burden to control the constraint not on the algorithm to be employed 
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but rather on the system of equations to be solved itself [TSfl^ . This is achieved by consider- 
ing an additional variable suitably coupled to the system through another equation so that, 
through the evolution, the constraint behavior is kept under control. While this method has, 
to date, received less attention than the constraint transport one, successful applications in 
diverse scenarios have already illustrated its usefulness (e.g. [T^fT^lHllTBlITT] ) . 

Regardless of the technique employed, boundary conditions can play a significant role and 
may spoil all efforts towards a correct implementation in the absence of boundaries. Clearly, 
even when a stable option is adopted, it need not be consistent with the goal of preserving 
the constraint (be it at round-off or truncation levels) and so carefully designed boundary 
conditions must be formulated. Commonly employed options include straightforward outflow- 
type conditions -which do not enforce the constraint- or "absorbing" boundary conditions 
which aim to reduce the influence of spurious effects induced at the boundary in the numerical 
solution [18]. 

In the present work we concentrate on formulating constraint preserving boundary condi- 
tions for the Newtonian ideal MHD equations (for systems with and without divergence clean- 
ing). Such a task is intimately related to the hyperbolic properties of the equations and so 
we re-analyze the system of equations and discuss alternatives for controlling the constraints 
through the divergence cleaning technique. 

We therefore organize our presentation along the following lines. In section II we analyze the 
system of equations and the formulation of boundary conditions beginning with a simplified 
model from which we draw the strategy to apply in the complete MHD system. Section III 
presents a series of tests that highlight the benefits gained by our construction. We conclude 
in section IV with some final comments. 



2. The Newtonian equations & constraint hyperbolic cleaning 

The equations describing the ideal Newtonian MHD equations in terms of the variables 
(p,p, v\B') are pni2j): 

dtp = -V.(pt;*), 

dtP = —vWiP — •ypVjV^ , 

dtB' = -Vj{v'B'-v'B'J). 

In this form, this system of equations is weakly hyperbolic since there is no complete set of 
eigenvectors. This indicates instabilities are likely to arise unless some modes are carefully 
controlled. The constraint transport technique attempts to do so at the algorithmical level by 
enforcing diB^ = at the discrete level. We here choose an alternative approach where the 
situation is remedied at the analytical level through the addition of an extra field coupled to 
the system in a suitable way P^fTl] . To simplify the discussion, we first consider a simpler 
model which captures essential features of the main system. 
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2.1. Preliminaries 



To illustrate this approach consider first the system obtained when the fluid field variables 
as given and stationary. The resulting system describes a simple model that shares key prob- 
lematic features of the original ideal MHD system ([T]), although in the latter case additional 
aspects enter under consideration. In this simple case one just has the magnetic field which 
evolves under, 

dtBi = -Vjiv^B' -v'B^) . (1) 

This equation is already weakly hyperbolic, since a plane-wave analysis indicates that when 
the wave number vector is perpendicular to the velocity vector there exists only two linearly 
independent eigenvectors. To see this consider the case where the wave number vector is normal 
to the velocity vector, and define coordinates axis such that the first axis lies along the wave 
vector -so ki = {k, 0, 0)- the second along the velocity vector, -so f * = (0, v, 0)- and the 
third perpendicular to both of them. Then the eigenvalue-eigenvector problem, obtained by 
expressing = jj^'^t+ki^'^ gives rise to the following problem: aU = MU (with a and U a 
frequency and a vector to be determined) and M defined as 



M 



/ o^ 

A;f 

\ oy 



(2) 



which is clearly non diagonalizable. Thus, there is no stability in the usual norm sense. 
Interestingly the divergence of the magnetic field is preserved by equation ([T]), that is, if one 
defines D := VkB^ and takes the divergence of ([T]) one obtains the induced evolution equation 
for D as dtD = 0. Hence, if D is zero initially it remains so as along as the time integration 
lines do not intersect a boundary. However, we note that if instead of considering the L2 norm, 
one considers a different norm defined by adding a term proportional to the divergence of Bi, 
namely a norm of the type: 



£:= J [B^ + c{VkB''Y]dV 



St 

one can show that S = and so this energy is controlled. Therefore, from an analytical point 
of view equation ([T]) is not a bad one. 

At the numerical level, things become more delicate as generic violations of this constraint 
will arise due to truncation or round-off errors which might grow unless a careful implemen- 
tation of the equations is adopted. Furthermore, even when an integration scheme is available 
that controls the constraint in the absence of boundaries, (which is not difficult to obtain in flat 
space in Cartesian coordinates -which we use in this work-), the boundary conditions modify 
the equations at the boundaries which can cause the constraint to grow and propagate to the 
interior. Thus, we next discuss a way to achieve a more robust behavior. Rather than doing so 
by the particular numerical algorithm employed we work first at the analytical level and adopt 
a strategy that would help even beyond the simple problem adopted here. 
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To remedy the lack of strong hyperbolicity and good constraint propagation a possible ap- 
proach can be adopted which makes use of the freedom to add to equation ([T]) any term pro- 
portional to that divergence. The proposed system is (in analogy with that in [14j): 

dtBi = -Vjiv^B' - v'B^) - ciVicj) , 

dt(p = -ciVkB'-s^. (3) 

Notice that for = the modified system is equivalent to the original one, and (p has 
as a source proportional to the constraint. Thus if boundary and initial data are such as to 
guarantee that V^-B^ = 0, and trivial data is given for (p the solution to the modified system 
is also a solution of the original equation. The advantage of this modification is that now the 
system is strongly hyperbolic, so it has a well posed initial problem and is stable irrespective 
of whether or not VkB^ vanishes. Furthermore, the induced evolution equation for D is no 
longer trivial, and it implies that the field D propagates with speed q. Consequently one can 
make the constraint violations to propagate away from the integration region, and, if correct 
boundary conditions are given, leave the computational region entirely. Additionally, the last 
term in the equation for in ([3]) induces a decay of the constraint for s > as it travels along 
the integration region further helping to keep it under control. 

Drawing from this exercise, other modifications are certainly possible; for instance, one could 
consider a more complex system given by, 

dtBi = -u^VjBi + B^VjUi - B.VjU^ + (1 - a)u'V jB^ - ciVi<p 
dtcj) = -(3u^Vj(j) - ciV'Bi - s(f) . (4) 

Which, for a range of values of the parameters a and f3 for is strongly hyperbolic. Among 
these, of particular interest is the system with a = (3 = 1, which is strongly hyperbolic even 
in the limit q — 0, so its associated initial value problem would be well posed even without 
the coupling to the field 0. Another interesting choice is the one a = f3 = since the resulting 
system can be expressed in conservative form. This property is often preferred in applications 
where shocks are present as one can easily take advantage of special algorithms defined to deal 
with these issues (see e.g. [T^2UII21f22ll23j ). 



Boundary conditions To fix ideas, we now study the possible boundary conditions for the 
above system (with a = (3 = 0). Our goal is to define conditions which can yield a well 
posed problem and, if possible, guarantee no violations are introduced into the computational 
domain. To this end we must determine which are the incoming and outgoing modes off the 
boundary surface as this information is key to understand what data is freely specifiable. For 
this we seek a solution of the form (5*, 0) = jj^'^t+^t^^ ^ where f/ is a vector to be determined 
along with the frequency a, and rii is the outgoing normal to the boundary under consideration. 
Applying this solution to the system we obtain: 

aB' = -B'vn + v'Bn - cin'(p , 
cr0 = , (5) 



4 



where Bn '■= B^rii, and f„ := v^rii. This problem of eigenvalues/eigenvectors can be solved to 
determine the incoming (cr > 0), tangential (a = 0) and outgoing (cr < 0) modes. A necessary 
condition for a well posed problem for hyperbolic systems indicates that data must be given 
only to the incoming modes, since the others are determined from the inside of the integration 
region [2321] • Recall that the incoming modes can be defined even as linear functions of the 
outgoing ones so long as the coefficients are small enough. 

At a given boundary, there are three eigenvalues, (ctq := —Vn, ct± := ±c;). The first corre- 
sponds to two linearly independent modes which span the tangent space of the boundary point 
(i.e. they are perpendicular to rii). They are positive whenever the velocity is incoming, these 
are the modes which are dragged along f The other eigenvalues correspond to the cases where 
we can choose i?„ := b arbitrarily and. 



B' 



(t>* ± n^ci)b 



Vn ± Cl 

Notice one can always adopt a value for c; larger than the maximum velocity expected so that 
denominator is never zero. The eigenvector corresponding to the positive eigenvalue determine 
the combination of variable that must be always suitably given to preserve the constraint. The 
complete expression for the eigenbase and its co-base is. 

f/°:=(0,ei,0) , a = -Vn, 
f/° := (0, 62, 0) , a = -Vn, 

Vn ± Cl 

where the generic vector is U := {Bn, B\ 0) with B'^rii = 0, and {ci, 62} are two orthonormal 
vectors on the tangent of the boundary. 
The co-base is given by 

VnVl CiVi , 



Cl -vl cf- vl 

a2._( ^nV2 C1V2 ^ 

cf - vl cf - vl 
^±:= ^(1,0,^1), 

where Vi := e\Vi, and f 2 := 63^4. In order to see how we must define boundary values consistent 
with the constraints we now analyze the subsidiary system determining the evolution of the 
constraint. In order to treat it as a first order system we define a new variable 5i := Vi0, whose 
evolution equation is determined by (the time derivative of) the evolution equation of 0. The 
full subsidiary system is then: 

dtD = -ciV'6, , 

dt6i = -ciWiD - s5i . (6) 

Here again we must investigate now which are the incoming and outgoing modes for this system. 
If we can impose boundary conditions to the main system so that we ensure no incoming mode 
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is created for this subsidiary system, then uniqueness of the (trivial) solution would guarantee 

that the solution of the main system has vanishing constraint quantities and so it is a solution of 
the original system. As done previously, plane wave solutions of the above system perpendicular 
to the boundary give rise to the following eigenvalue/eigenvector problem: 

kD = —cin'Si , 

K6i = -ciTiiD , (7) 
and its eigenvalues are {kq — 0,k± — ±q} with the corresponding eigenvectors: 





:=(0,ei,0) , 


= 0, 




(0,62,0) , 


«o = 0, 




■={h- ^ 


, =Fl) , K± = ±Ct 



where the generic vector is V := (5„, 5\ D) with 5'"ni — 0, and {ei, 62} are two orthonormal 
vectors tangent to the boundary. 
The co-base is given by 

0^=(O,ei,O), 

e^= (0,62,0), 

e±:= ^(1,0,^1). 

Here again the theory of boundary conditions indicate we should give as a boundary condition 
the incoming mode as a linear function of the outgoing one. That is we should ensure that the 
boundary data for the main system is such that: 

e+(y) + ae_(y)=o , |a| < i (8) 

This is a linear combination of space derivatives of the fields and can be implemented in several 
different ways depending on the problem at hand. Here we choose to implement it at the level 
of the evolution equations, namely by modifying the evolution equations at the boundary so 
that the constraint is satisfied there. From the evolution equations of the main system we can 
solve for D and 5„ in terms of the time derivatives of S„ and (p (and the remaining equations). 
Hence 

5n = {dtBn - Fn) , 

Q 

D^—{dt<j>-F^), (9) 
Q 

with Fn = —nidj{u^ B"^ — u^B^), and F^ = — s0. 

We can then translate the above boundary condition into: 

- Bn{l + a) + 0(1 - a) + Fn{l + a) - F^{1 - a) = , (10) 
which can be re-expressed as: 

- 2 (e+(U) + ae^U)) + F„(l + a) - F^(l - a) = . (11) 
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For the particular case of setting the incoming constraint modes to zero (a = 0) this imphes, 
0AU) = liF^-F^) . (12) 

Thus, the boundary condition defined by the equation above is not only maximally dissipative 
but also enforces the constraint. 

2.2. The complete system 

We now turn our attention to the complete system and develop an analogous strategy fol- 
lowing the discussion considered above. In most cases the literature on the subject considers 
the MHD system expressed in terms variables U = [p, e, m\ B^) rather than U = {p, p, u^, B^) 
-ie the internal energy e instead of the pressure p-. While this involves a simple change of vari- 
ables it turns out the characteristic decomposition is simpler in the latter case. So, our evolu- 
tion system will be given by the former, and we will obtain the characteristic structure for the 
latter followed by a straightforward change of variables. Our system of interest is a variation 
of the initial MHD equations as (see 



dtp = -V^ipv'), 
pdtv' = -pv^Vjv' - V'p - BkiV'B^ - V'^B') - aB'WkB^ , 
diB' = -Vj{u^B' - u'B^) - au'VjB^ - ciV'cj) , 

dte = -Vi((e + P + ^B^y - B'v ■ B) - av'BiWuB'' - ciB^Wkct^ , 

dt(() = -au^Vj(f) - ciVjB^ - s(j) , (13) 
where 

p := (7 - 1) (e - Ip.^ - ^5^) , c^= ^ (14) 

Thys system includes both the possibility of divergence cleaning (q 7^ 0) and Galilean invari- 
ance {a = 1) and is strongly hyperbolic, hence has a complete set of eigenvectors. This will 
allow us to introduce a new boundary treatment which ensures no constraint violations are 
introduced through the computational domain boundaries. The parameter a controls the free- 
dom of adding the constraint equation. In what follows, we shall use the values a = 0, 1. The 
value a = corresponds to the purely conservative system (notice that velocity equation can 
be re- written with the help of the first equation as momentum conservation). The case with 
a = 1 gives rise to a system that is Galilean invariant (see [H]) and also to the case where the 
system is strongly hyperbolic [26|27J irrespective of the field 0. Hence one can study the limit 
where the (j) field decouples and the constraint propagates along the velocity field. 



Boundary system As done previously, we now consider a boundary with outgoing unit normal 
n* and obtain the characteristic decomposition at this boundary. For notational purposes we 
will employ an overbar to denote the perturbation on a given variable, for instance p will 
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denote the perturbation of p which will be considered as a fixed background quantity. The 
characteristic decomposition is then determined from the system, 



Crp = -Vnp- pVn , 

avi = -VnVi + BiBn/p+ {Bn/p)B, - ni{BjB^ +P)/p, 

aBi = -VnBi + ViBn - BiVn + BnVi + UiCj) , 
(^P = -clpVn - (7 - l)BjV^Bn - Vnp + v{l - '^)Bn(t) ■ 

CT^ = c^Bn (15) 

The solution to this system is cumbersome and requires dealing with different cases. The full 
analysis and solution is presented in appendix El and expressed in terms of a vector whose en- 
triesare:[/ = (p, Vn, Vj, B^, Bi, p, 0)-^indicatingperturbationsoff (p, f^rij, _B*nj, 0)-^ 
respectively with v\ B^ indicating the components of f * and -B* orthogonal to the boundary. 
While the full characteristic decomposition is lengthy, the characteristic decomposition asso- 
ciated to Bn and is particularly simple. In fact, the associated co-basis is given by 

0L± = 7T^(O,O,O,l,O,O,±cr'), (16) 

i.e. involving essentially just Bn and 0, hence the imposition of suitable boundary conditions 
will be directly related to that of our discussion of the simplified model in section [2711 



2.3. Boundary conditions 



Having the complete characteristic structure it is now rather straightforward to determine 
the possible boundary conditions. A simple one is to set all incoming modes to zero but this, 
generically, will not be consistent with the constraint. To this end, as explained in section 
12. H one must analyze the induced evolution for the constraint and obtain from it a recipe for 
what to provide to the incoming modes. In the general case discussed here, notice that Ql± is 
essentially the same that gives rise to equation (JT2l) since is non-trivial only for the Bn and A 
components. As a result, one can follow the same strategy and formulate constraint preserving 
boundary conditions by enforcing 

if = UL+ - F„)) = £iU) , (17) 

where £{U) denotes the (maximally dissipative) constraint preserving boundary condition 
defined by equation f|T7j) . 



3. Numerical Tests 



In this section we present the results of tests aimed to examine the solution's behavior with 
the different possible choices. To this end we constructed an implementation of the MHD 
equations (eqns. [T3|) on a 2-dimensional domain with Cartesian coordinates (x, y) G [—L, L] x 
[— L, L]. To discretize the system we employ a set of techniques which guarantee the stabihty of 
generic linear hyperbolic systems. These techniques are constructed to satisfy at the discrete 
level, the conditions holding at the analytical one to prove well posedness of a problem through 
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energy estimates [23]. Thus we adopt discrete derivative operators satisfying summation by 
parts, a 3rd order Runge-Kutta time integration and add dissipation through a Kreiss-Ohger 
type operator consistent with summation by parts. Finally, in all the tests considered we add 
a homogeneous 'atmosphere' throughout the computational domain in the initial data (and a 
'floor' during the evolution) for the density and the pressure to prevent negative values from 
arising in the simulation. We typically adjust the atmosphere (floor) values to be 10~^ (10~^) 
times the maximum of the initial values of p and e respectively. Notice that this straightforward 
implementation is not expected to handle shocks or discontinuities, however, since our goal is 
to examine boundary conditions we restrict here to scenarios where these effects do not arise 
within the time of interest. 

For testing purposes we adopt a few different initial scenarios which are based on modifica- 
tions of two commonly employed initial data sets: 



Rotor test: This is a variation of the MHD Rotor problem [28] which is smooth and allows 
for considering an additional flow field and initial constraint violation. 



10 if r < ro 

1 + 9/(r) if ro < r < ri 

1 if ri < r 

y 



(18) 



w''(x, y) = i/^' + < 



^^(x, y) = z/^ + < 



and 



-f{.'r')vo- 
ro 

-f{r)vo- 
r 





f{r)vo— 
rp 

fir)vo- 
r 







if r < To 
if ro < r < ri 
if ri < r 

if r < ro 

if ro < r < ri 

if ri < r 



(19) 



(20) 



p = l 

3^ = 
d) = 



(21) 
(22) 
(23) 
(24) 



with ro = 1, ri = 2, r = y/x^+y^, 7 = 1.4, vo a constant and /(r) an interpolating function 
For the general tests we adopt 

(ri — r) 



fir 



(25) 



(ri - ro) ' 

However, for the convergence test we adopt a smooth interpolating polynomial and velocity 
profiles as 
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/(r)=P5(r) (26) 
Vo = vy/{l + rY (27) 

with P^{r) a polynomial of order 5 such that P^{rQ) = 1; P^iri) = and both first and 
second derivatives at r = {tq, ri} are zero. Thus, the corresponding initial data is at least 
throughout the computational domain. 

Finally, the functions {u^, v^) and are included to consider further modifications for specific 
tests, in particular choosing : 

- = Ke~^ allows one to introduce data violating the constraints. 

- (i/^, z/^) = — e (l — e"'' ) (x, y) defines data infiowing to further test boundary issues. 



Blast test: This is a variation of the MHD Blast problem [29] . 
The initial data is given by. 



p 


= 1.0 




= 0.0 


Vy 


= 0.0 


Bx 


= 4.0 


By 


= 0.0 


P 


= e-^' 




= 0.0 



(28) 

Additionally, we also studied the solution's behavior employing slight variations of this data set. 
In particular, cases with initial incoming or outgoing velocity fields at points of the boundary. 
In this way we could test different modes on the boundary which become in some cases outgoing 
or incoming according to the velocities chosen. 



Boundary conditions adopted For boundary conditions we adopt one of the following three 
possible cases: 

(i) Freezing boundary conditions. Defined by Ui = 0, for i = 1..7 (denoted by FR). 

(ii) Outfiow boundary condition (incoming modes set to zero). Defined as S+f/ = 0, with 
S+ = Ei Uj'Oj, for j such that > (denoted by NI). 

(iii) Constraint preserving boundary condition. Defined as S+f/ = £{U), as defined by eqn. 
([HD (denoted by CP). 

in all tests, and compare the solutions' behavior when employing each option. We begin by 
considering the fiux-conservative form of the equations (keeping a = 0) and then examine 
particular cases with the Galilean invariant form {a = 1). 

3.1. Testing the implementation 

In the first test we confirm the overall convergent behavior of the numerical solution when 
considering sufficiently smooth initial data. We evolve the version of the Rotor initial data 
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Fig. 1. Behavior of the L2 norm of the difference between solutions obtained with gradually improved resolutions. The difference 
between them converges to zero as expected. 

(with no constraint violation or background fluid flow) for three different resolutions A; = 
1.5/2' (/ = 0, 1, 2) and check the pair-wise difference of the numerical solutions obtained de- 
creases as expected. Figure [1] illustrates the behavior of ||F(Ai) — F(A/_|_i)||2 (with F either 
p or B^). As is evident in the figure, as resolution is improved the differences decrease as ex- 
pected. For all remaining tests we present results for the finest grid employed with 1 = 2. 

3.2. Blast initial data 

In this test, we adopt the blast initial data, evolve it for different choices of boundary condi- 
tions setting a = 0, Q = 20 and s = 1 and examine the constraint's behavior in each case. Fig- 
ure [2] shows the L2 norm of the constraint for the different boundary value conditions. Clearly 
the numerical solution obtained with constraint preserving boundary conditions is superior 
by about an order of magnitude in constraint violation than the no-incoming case and almost 
three orders better than that obtained with the freezing boundary condition. An important 
point to emphasize is that a closer inspection of the solutions obtained reveals that the main 
contribution to the error originates at the boundaries in all cases (though with essentially the 
same behavior as far as the error's magnitude with respect to the boundary condition adopted) . 
The norm displayed in Figure [2] is calculated over the whole computational domain ignoring 
the last two points at all boundaries to avoid placing excessive weight on the violation at the 
boundary. Nevertheless, as indicated, constraint preserving boundary conditions give rise to a 
solution whose constraint is violated the least. 

3.3. Rotor initial data. Effects of boundary conditions and divergence cleaning 

We adopt this data to further examine the solution's behavior and allowing for initial viola- 
tions of the constraint. To this end we adopt /t = 10~^. 
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Fig. 2. Behavior of the L2 norm of the constraint for the different boundary conditions. The violations in each case grow from 
round-off values before settling to an approximately constant value. The values throughout the run for this norm improves as 
the boundary condition adopted is refined as expected, in fact the solution obtained with the constraint preserving boundary 
condition preserves the constraint by at least three-order of magnitude better than the one obtained with the freezing condition. 

Figure [3] illustrates the solution's behavior with the different boundary conditions and with 
and without the use of the divergence driver. The addition of the divergence driver allows for 
a dynamical reduction in the constraint violation until t ~ 4.6, at this time the propagating 
modes interact with the boundaries which become the main source of error. Here again one 
sees that the constraint preserving boundary condition gives rise to significantly smaller er- 
rors in the solution. Interestingly however, adopting the no-incoming modes together with the 
constraint damping field provides a reasonably similar behavior. This is due to the fact that 
the no-incoming boundary condition allows the outgoing constraint violating mode to leave 
the computational domain, while the incoming constraint violating mode, generated by not 
imposing the constraint preserving boundary condition, is damped to a significant degree by 
the divergence cleaning technique. 

Another interesting behavior is revealed when varying q. This affects the constraint violating 
modes' propagation speeds which in turn has a strong influence in the solution's constraint 
behavior. In what follows we adopt the constraint preserving boundary condition. For this 
test we use the Galilean invariant system, setting s = (so as not to include damping) and 
choose Q ranging from 10 to 800 together with adjusting the time step accordingly in order 
not to violate the Courant condition. Figures [Hand [5] (with initial constraint violation) show 
the effect of varying q in two cases, the Rotor initial data and its modification to include a 
background velocity given by (z/^' = ztO.lx, u'^ = ±0.1y). We concentrate on the latter case 
for there the effect is more striking. Notice that, in figure HI the initial plateau corresponds 
to round-off values since the constraint is initially satisfied to that level. Once the non-trivial 
part of the solution reaches the boundary a significant violation of the constraint is generated. 



* For smaller values of cj (c; < 10) the code became unstable due to instabilities generated at the corners of the computational 
domain. 
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le+00 




Fig. 3. Behavior of the L2 norm of the constraint for different options with some non-trivial initial violation of the constraints. 
As is evident in the figure, the divergence cleaning is able to damp the constraints by several orders of magnitude while the 
violation is present in the bulk of the computational domain. After the solution reaches the boundary, the boundary values 
induced there dominate the violation of the constraint and again the behavior is significantly improved by the no-incoming 
and constraint preserving conditions. 




Fig. 4. Behavior of the L2 norm of the constraint for different values of the coupling constant c; with constraint preserving 
boundary condition. The norm is calculated over the whole computational domain (whole) or over an interior region of 3/4 its 
size (masked). The latter illustrates how, despite not including a region close to the boundaries, boundary induced effects are 
evident throughout the domain. As the value of c; is increased a slight improvement is achieved. 

Depending on the boundary condition, and on the characteristic velocities of the system under 
consideration that violation propagates to the inside or just leave the domain. In our case, the 
ability of the constraint preserving boundary condition to allow constraint violating modes 
to leave the domain, results in smaller errors for faster propagating cases. This is a result of 
violations induced by the evolution leaved the computational domain more rapidly. 
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Fig. 5. Similar to figure |4] but with initial data containing constraint violations. 



4. Conclusions 



We have investigated the numerical solution to the (ideal) Newtonian magnetohydrodynam- 
ics equations and developed boundary conditions which at the continuum level are constraint 
preserving and so would not introduce further violations in the computational domain. At 
the numerical level these conditions do introduce some truncation-error level violations which 
converge away with resolution. These boundary conditions developed are consistent with max- 
imally dissipative ones and so the discrete energy of the system remains bounded. We also 
examined the constraint's behavior when enlarging the system so as to couple an extra field. 
The addition of such field, together with a suitable modification of the equations induces a 
non-zero speed in the propagation of the constraint and allows for driving its value to zero. 
We have studied, in particular, two of the many equivalent systems. One which is fully conser- 
vative -which is only weakly hyperbolic without the addition of the extra field which renders 
it symmetric hyperbolic- and a Galilean invariant one which is symmetric hyperbolic -even 
without adding the extra field- but is not expressible in conservative form. 

In order to examine the solution's behavior we implemented the equations in Cartesian 
coordinates and with a spatial discretization that preserves the initial constraint error. This 
allows us to separate bulk from boundary effects and observe that the main violations do indeed 
occur at boundaries. To examine the boundary condition effects on the evolution we adopt 
three types of boundary conditions: 

- The most direct one, and crudest, sets to zero all right hand sides in a buffer boundary region. 
This might lead to inconsistencies for one might be prescribing conditions to outgoing modes. 
This condition is essentially equivalent the most commonly employed one of flux-copying 
near boundary points. 

- The second one is the "no-incoming" boundary condition where one projects out all incoming 
modes in the evolution equations leaving the rest (tangential and outgoing modes) intact. 
The implementation of this condition is slightly delicate as modes can change directions over 
time, and so may turn from being incoming to being outgoing and vice- versa. 
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- The third condition preserves the constraint in the sense that is deduced by setting to zero 
the incoming mode related to the constraint propagation. This condition is a modification 
to the previous one where now instead of projecting out all incoming modes, one incoming 
mode is fixed so that the incoming constraint mode is zero as a result. 
Not surprisingly the first option displays the worst behavior for the constraint as not only 
nothing is done to minimize the constraint violations introduced through the boundaries but 
further inconsistencies in the solution are induced there. As a result, in the best case the 
constraint mode bounces back from the boundary into the integration region. The second 
option performed substantially better displaying a gain of an order of magnitude in constraint 
preservation. This is a result of the boundary condition allowing for one of the modes - 
carrying constraint violations- to leave the integration region. With the third alternative an 
additional order of magnitude (at least) is gained with respect to the non-incoming condition. 
This condition not only allows for constraint violations to leave the computational domain but 
does not introduce significant violations at the boundary. 

Thus we see that appropriately handling the boundary conditions the problem of constraint 
behavior can be significantly controlled. It should be mentioned that we still do not have 
a mathematically rigorous proof that the constraint preserving boundary conditions is well 
posed. However, in simpler systems -like the toy model we discussed at the beginning of the 
paper- a proof could be devised at least for the case where the eigenvalues do not change sign. 
Further work in this direction is needed to put the analysis and results obtained in stronger 
grounds. 
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Appendix A. Characteristic decomposition 

In this appendix we revisit the characteristic analysis of the ideal MHD system coupled to 
the divergence cleaning field (a related discussion in the absence of this field can be found in 
[30;31j). Using the decomposition on normal and tangential parts, and defining d := cr + we 
get: 

dp = -pvn , 

dvn = 2BjpBn - l/pB,B^ - 1/pp, 

dBn = V„Bn + , 

dp = -clpVn - (7 - l)BjV^Bn + (1 - l)Bn(j) , 
dBv = B^/pBn-Bn/pP, 

dBB = BvBn - B^Vn + BnBv + Bn(f) , 
di)i = Bi/pBn + Bn/pBi, 
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dBi = ViBn - BiVn + BnVi . 

<t4> = cfS„ (A.l) 

For the sake of clarity we will present the solution to the eigenvalue/eigenvector problem 
along different cases which are naturally divided by the behavior of certain modes. 



A.0.1. BASIS 

• Normal modes. These correspond to ^„ = = 0, where the equations become: 
dvn = -BB/p - p/p , 

dp^-clpVn, 

dBv = -{Bn/p)p , 
dBB = -B^Vn + Br^Bv , 

d^i={Bn/p)Bi, 

dBi = -BiVri + BJi . (A.2) 

Setting d = one can show that all but the first component must vanish while itself can have 
any value. So a first eigenvector is given by 

C/o = (l, 0, 0, 0, 0, 0, 0), (A.3) 

where the entries are: U = {p, Vn, Vi, Bn, Bi, p, 4>). 
Using all the equations we get the eigenvalue condition: 

d'-{cl + Byp)d' + clBl/p = 0. (A.4) 

Prom which we get four solutions: 



dp± = ±y {{cl + Byp) + ^{cl + B^pY - A{Bl/p)ci) /2 (A.5) 
dm± = ±^ ((c2 + B^p) - ^{cl + B^pY - A{Bl/p)4) /2 (A.6) 

Notice that when ^ the dp solutions tend to ic^, while the dj„ — > as ||i?„||/y^. The 
dp solutions tend to the fluid solutions while the d^, to pure magnetic ones, so they decouple 
in this limit. One thus must choose the eigenvectors' normalization in a suitable way to reflect 
this behavior. For the dp solutions we set = 1, and use the last two equations to compute 
the normal magnetic fleld and the normal velocity, obtaining: 

UP^ = U-^, 1, 0, T^Bi, t'^, o] , (A.7) 

where Rp :— dp — B^/ p, has a flnite hmit as — > 0. 
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For the dm solutions we proceed the other way around, we fix BB = BiB = BiM\ where 
Mi = Bi/\B\ and so compute Vn, and p from the second and third equations of system flA.2p . 
obtaining: 

fBB TdrnBB ±Bn c'iBB \ , 

f/M± = ^ , —J^M,, 0, M„ ^— , , (A.8) 

where Rm '■= df^ — has a finite hmit as B"^ 0. 

In two dimensions these are all eigenvalues-eigenvectors pairs in the normal sector. 

In three dimensions we can choose Bi = Mi perpendicular to Bi, but otherwise arbitrary, 

that is BB = BiB = 0. In this case all other components vanish except the tangential velocity 
which can be expressed as Vi = da/BnAi with A = {—B2, Bi)/\B\ and da± = -^Bnj ^ (so 
that the last two equations can have a nontrivial solution). The corresponding eigenvector is: 

UA^ = (^0, 0, ±^A, 0, Ai, 0, oj . (A.9) 
• (j) modes 

We now look at the modes which are induced by the introduction of the field (p. There, from 
the third and last equation of flA.ip we get a 2x2 system which implies: a = ±q and so, d := 



di = ±ci — Vn- If we fix Bn = 6 we can obtain the scalar components by solving the following 
subsystem of ( lA.ll) : 

diVn + BB/p + p/p = 2(5„/p)6 , 

clpvn + pdi = (-(7 - l)Bv + (1 - 7)a5„)6 , 
-B^v„ - diBB + BnBv = -{Bncr + Bv)h , 

{Bjp)p + diBv = {B^/p)h. (A.IO) 

From this subsystem we get: 

Vn = [-{d^i - BIIp)F2 + d}F; + dip{d}Fi - {Bnlp)F^)]/6i , (A.ll) 
p={-clpVn + F2)/di, (A.12) 

where: 

5i = p{dt-d]{cl + B^/p) + clBl/p), 
-2BJ) 

P 

F2 = {l--f){Bv + Bna)b, 

F3 = -F2/(l-7), 

B% 

Fa = . 

P 

Once we have Vn we can solve for the vectorial components of the system: 
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di^i-{Bjp)Bi = Bib/p, (A. 13) 

diBi - BnVi = Vib - BiVn , (A. 14) 

obtaining: 

= [vibBjp + Biidib - v„B^)/p]/S, , (A. 15) 

Bi = [vidib + Bi{-diVn , +bBl/p)]/Ss, (A.16) 
where Sg '■— df — B^/p. Combining the above intermediate steps we finally obtain, 

C/L± = (^-(p/di)vn, [-{df - Bl/p)F2 + d^'Fa + d^F^ - ^)]/Si, 

bB^ 

[VibBn/ p + Bi{dib - VnBn)/ p]/Ss, b, [vidib + Bi{-diVn ^ -)]/^s, 

P 

t^^^. (A.17) 



A.0.2. CO-BASIS 

Net we compute the co-basis which we will use to construct the suitable projector for en- 
forcing different boundary conditions. As in the previous analysis, we divide our task along 
different eigenvalues for simplicity 

• eo 

Since Uo has only one non-vanishing component (the first one) all elements of the co-base, 
except the corresponding one have a vanishing first element. 
The first element, is simply: 

Go = (1, 0, 0, A, 0, -cf, B). (A.18) 

To find the remaining elements we notice that if we define: VP := {U L+ + U L_) and V M : = 
{UL^ — U the first has a zero in the last component (corresponding to 0) while the second 
has a zero in the fourth component (corresponding to Bn) ■ Thus contraction of ©o with VP 
and VM will leave either A or B only as unknowns. Thus we might proceed to setting A and 
B temporarily to zero, perform the contraction and extract what these must be obtaining, 

A^^, (A.19) 

B^-^f^\ (A.20) 
2bci ^ ' 

where we first temporarily set to zero A, and B uiQq. 

•©p± 

Prom 
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ep±(C/M+ - C/M_) =0 and 
ep±(C/M+ + C/M_)=0 

we get that the following structure for it: 

Qp± - U C, C^, D^, E^, . (A.21) 
Now, from 1 = Q p±{U Pj^ + U P_) we obtain: 

C = ^p3ii1^ (The same for both) . (A.22) 

2{RpRm-dlB^/p) 

While from ±1 = 0p±(t/P+ - UPJ) we deduce 

^ ±RpR„.dj, ^ 

2cl{dlB^-pRpRm) 

The remaining two components can be computed using VP and VM as above obtaining, 

(A.24) 

_-0p±(VM) 

2b^, ' ^^-^^^ 

where we set to zero D, and F in 0p±, as in the previous case. 

• ©Mi 

From Qm±(UP+ - UP^) = and eM±{UP+ + UP J) = we get that the following structure 
for it: 

©Mi = fo, ^^4^, iVi, 5„ L^ t\ . (A.26) 

\ clRpp j 

Prom 1 = ©M±(C/M+ + C/M_) we obtain: 

U = , ^^-^p^^ . (A.27) 

2iB^dl/p-R^Rp) 

From ±1 = 0M±(f/M+ - [/M_) we obtain: 

Ni^^ ^^-t^-^^p^' . (A.28) 

2Bjp{B^dlJp-R„,Rp) 

The remaining two components, S and T arc computed as in the previous case, using: VP 
and VM removing the 5„ and the last component (f) in 0m±- We get, 

S,^--^^, (A.29) 
r.^=%^. (A.30) 

20Q 
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• eA± 

We get the co- vector of [7^4 by direct computation: 

Qa± = {o, 0, G±, 0, , (A.31) 

where the components, G and if are given by: 

G,.Z^, (A.32) 
removing the components Bn and in Ga±- 

It remains now to determine the last two elements. Since the first 7 eigenvectors span completely 
the seven dimensional space given by (S„ = = 0) the co- vectors have only components in 
that subspace and are given by: 



Caveat: Singular points 

The characteristic structure outlined above may change as some eigenvalues can change mul- 
tiplicity for particular values of the fields. These special cases require further analysis. Recall 
the eigenvalues are given by: 





= -v„ 




a+ 






ap 


dp 


-Vn: 




dm ~ 


■ Vn , 




dm 


- Vn 




= da- 


Vn, 




= -da 


-Vn- 




= Cl, 






= -ci, 





they can cross when dp = 0, w„ = ±q, dm = Q and dp = dm- The first two cases can be dealt 
with by appropriately chosen parameters so that they would not occur in physically relevant 
scenarios. For instance, the first case would imply a vanishing sound speed and will not be 
considered since we will not deal with a fliiid describing dust. The second case could be avoided 
by simply taking q large enough. Therefore, we concentrate on the other two: 
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For this we need that 





4c^^ (A.35) 
P 



dl = -t^,- (A.36) 



which only happens when Bn = 0. Near that point we have, 

pel + 52 

At this point we compute exphcitly the eigenspace for the coinciding eigenvalues and choose 

eigenvectors with good limit. We use them in a sufficiently small neighborhood of this point 

using conditional statements in the code at every point in the boundary. 

In three dimensions, in this case, also da = 0, therefore the eigenvectors UM± and UA± are 

degenerates. 

* dm dp 

For this to happen we need 

= « - = [cl + ^ j - AclBl = [cl -^j + 2cA + ^. (A.37) 

Thus. B'^/ p = cj., and Bi = 0. Notice that in that case we have Rp = Rm = 0. We proceed in a 
similar fashion as above. Notice that when eigenvalues coincide all what enters in the boundary 
condition is the projector on that subspace, so one just has to choose the eigenvectors so that 
numerically that projector is robust in a small neighborhood. 

In the case of three dimensions, = B\l p. Hence, da = dm = dp and the eigenvectors UP±, 
U M± and U A± become degenerates. 
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